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ABSTRACT 


The  system  of  partial  differential  equations  governing  the  nonlinear 
transient  motion  of  a  cable  immersed  in  a  fluid  is  solved  by  finite  difference 
methods.  This  problem  may  be  considered  a  generalization  of  the  classical 
vibrating  string  problem  in  the  following  respects:  a)  the  motion  is  two 
dimensional,  b)  large  displacements  are  permitted,  c)  forces  due  to  the 
weight  of  the  cable,  buoyancy,  virtual  inertia  of  the  medium  and  damping 
or  drag  are  included,  and  d)  the  cable  is  assumed  to  be  nonuniform.  The 
numerical  solution  of  this  system  of  equations  presented  a  number  of  interesting 
mathematical  problems  related  to:  a)  the  nonlinear  nature  of  the  equations, 
b)  the  determination  of  a  stable  numerical  procedure,  and  c)  the  determi¬ 
nation  of  an  effective  computational  method.  The  computation  is  programmed 
for  a  high-speed  calculator  (UNIVAC  system).  The  solution  of  this  problem 
is  of  practical  significance  in  the  calculation  of  the  transient  forces  acting 
on  mooring  lines  due  to  the  bobbing  up  and  down  of  ships  during  the  period 
preceding  large  scale  explosion  tests,  as  well  as  in  many  other  applications 
involving  mooring  or'  towing  operations. 


1.  INTRODUCTION 


This  problem  arose  as  a  result  of  an  urgent  requirement  by  the  Navy 
in  connection  with  a  series  of  nuclear  explosion  tests  which  were  conducted 
in  the  Pacific.  In  preparation  for  these  tests  a  number  of  ships  were 
instrumented  and  moored  at  specified  locations  from  the  explosion  point. 
These  positions  had  to  be  maintained  intact  during  the  period  preceding  the 
explosion.  However,  the  bobbing  up  and  down  of  the  ships  due  to  ocean  waves 
during  this  period  could  excite  sizeable  transient  forces  in  the  mooring  lines 
which  might  break  these  lines  and  thus  result  in  the  loss  of  information  from 
the  tests.  Several  months  prior  to  these  tests  a  request  was  made  to  the 
Applied  Mathematics  Laboratory  to  calculate  the  magnitude  of  the  forces 
acting  on  the  mooring  lines,  for  waves  of  varying  amplitude  and  frequency. 

It  is  gratifying  to  report  that  in  spite  of  the  theoretical  complexities  of  this 
problem  and  the  absence  of  any  known  solutions,  the  Applied  Mathematics 
Laboratory  was  able  to  obtain  the  required  results  in  time  for  use  during 
the  scheduled  tests.  The  two  factors  which  made  a  theoretical  solution 
feasible  at  this  time,  whereas  it  would  not  have  been  possible  several 
years  ago,  were:  a)  the  availability  of  a  high-speed  computer  and  b)  the 
recent  progress  made  in  the  understanding  and  development  of  numerical 
methods  for  the  solution  of  systems  of  partial  differential  equations. 

Whereas  the  solution  to  this  problem  was  carried  out  as  a  result  of 

one  specific  requirement,  it  is  more  useful  to  regard  it  as  the  general 

problem  of  the  two-dimensional  motion  of  a  cable  or  rope  immersed  in  a 
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fluid.  From  this  point  of  view  it  may  be  considered  as  a  rather  broad 
generalization  of  the  classical  vibrating  string  problem,  and  it  becomes 
immediately  apparent  that  its  solution  is  applicable  to  a  wide  class  of 
engineering  problems  involving  the  motion  of  cables,  such  as  a)  the 
laying  of  intercontinental  telegraph  cables,  b)  the  towing  of  a  ship  or 
other  object  in  water,  or  c)  the  snapping  of  telephone  wires  as  a  result 
of  transient  forces  caused  by  storm.  This  problem  may  be  stated  abstractly 
as  follows:  Given  the  initial  conditions  (i.  e. ,  position  and  velocity  at  any 
time,  to)  and  boundary  conditions  (positions  of  end  points  at  all  times)  of 
a  cable  immersed  in  a  fluid,  determine  its  subsequent  motions.  The  motions 
are  assumed  to  take  place  in  two  dimensions. 

Forces  that  .are  assumed  acting  on  the  cable  are:  a)  forced  motion  of 
the  extremities  (end  points)  of  the  cable,  b)  damping  or  drag  as  it  moves 
through  the  fluid,  c)  virtual  inertia  due  to  the  motion  imparted  to  the  fluid 
d)  weight  of  the  cable,  and  e)  buoyancy.  Variations  in  the  mass  as  well  as 
other  physical  properties  of  the  cable  along  its  length  are  allowed.  The 
displacements  may  be  large  and  the  motions  rapid.  In  the  present  solution 
it  is  assumed  that  the  cable  is  inextensible  (cannot  be  stretched).  In 
subsequent  work  the  authors  have  carried  out  solutions  for  cables  with 
elastic  properties.  The  motions  are  not  restricted  in  any  manner  (except 
that  these  take  place  in  two  dimensions). 

The  solution  was  carried  out  by  the  method  of  finite  differences. 

This  method  consists  simply  in  replacing  the  derivatives  of  various  order 
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in  the  differential  system  of  equations  by  equivalent  ratios  of  finite  increments. 

This  substitution  results  in  a  system  of  difference  equations,  which  are 

algebraic  in  form,  and  hence  more  easily  tractable.  However,  in  order  to 

represent  a  valid  solution,  the  system  of  finite  difference  equations  so  derived 

must  possess  certain  mathematical  stability  (and  convergence)  properties. 

It  must  have  the  property  that  its  solution  progressively  increases  in  accuracy 

as  the  size  of  the  time  increment  used  in  the  above  representation  is  gradually 

decreased.  Unfortunately,  the  system  of  finite  difference  equations  initially 

proposed  for  the  solution  of  this  problem  did  not  satisfy  these  stability 

requirements.  A  search  for  a  stable  finite  difference  system  as  well  as  for 

an  effective  method  for  solving  the  resulting  system  of  finite  difference 

equations,  which  was  nonlinear  in  character,  added  to  the  complexity  of 

* 

the  problem. 
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2.  DERIVATION  OF  EQUATIONS  OF  MOTION 

The  problem  under  consideration  is  a  generalization  of  the  classical 
problem  of  the  motion  of  a  vibrating  string.  In  at  least  one  respect  the 
formulation  of  the  problem  given  here  will  not  differ  from  that  of  the  original. 
Specifically,,  we  wish  to  deduce  the  approximate  motion  of  a  steel  cable 
without  having  to  involve  ourselves  in  the  explicit  computation  of  the  elastic 
forces  which  act  on  the  cable.  However,  the  formulation  will  depart  from 
the  original  in  a  number  of  respects,  namely: 

a)  Longitudinal  as  well  as  transverse  motions  of  the  line  must  be 
taken  into  consideration. 

b)  The  occurrence  of  large  displacements  from  the  equilibrium 
configuration  of  the  line  must  be  permitted. 

c)  The  weight  of  the  cable  must  be  taken  into  account  because  the  line 
may  stretch  from  one  level  to  another.  Thus,  even  when  the  line  is  in  static 
equilibrium,  the  tension  will  not  be  uniform  nor  will  the  line  be  straight. 

d)  Since  the  cable  is  submerged,  the  static  forces  must  include  the 
buoyancy  of  the  medium  and  the  dynamic  forces  must  allow  for  the  virtual 
inertia  of  the  medium.  Furthermore,  it  is  desired  to  make  provision  for 
damping  forces  due  to  the  drag  on  the  line  whenever  lateral  motion  is 
occurring. 

e)  Finally,  it  is  desired  to  suspend  concentrated  loads  at  one  or  more 
points  along  the  line  and  to  change  the  linear  density  of  the  cable  at  specified 
points. 
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The  best  approach  to  the  solution  of  a  problem  with  such  general 
specifications  appears  to  be  a  numerical  method  based  on  finite  difference 
approximations.  Inasmuch  as  we  are  committing  ourselves  to  the  eventual 
use  of  differences  in  both  the  time  and  space  dimensions,  it  will  be  simpler 
to  introduce  the  spacewise  discreteness  into  the  original  formulation  of  the 
problem.  We  therefore  proceed  at  once  to  the  derivation  of  the  equation  of 
motion  of  a  simplified  model  in  which  the  distributed  mass  of  the  cable  has 
been  replaced  by  a  series  of  discrete  masses  mj  attached  to  a  weightless, 
inextensible  line.  This  leads  to  a  system  of  ordinary  differential  equations. 
It  may  be  shown  that,,  in  the  limit,  the  resulting  equations  pass  over  into 
the  corresponding  partial  differential  equations  for  the  motion  of  a  sub¬ 
merged  cable. 
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Figure  1  -  Mooring  line 


Figure  2  -  Discrete  representation 
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Figure  1  shows  a  typical  configuration  of  the  system  with  the  cable 
attached  to  a  float  at  the  surface  and  anchored  to  the  bottom.  Also,  a 
heavy  load  is  suspended  from  a  point  near  one  end  of  the  cable.  Other 
boundary  conditions  are  possible,  but  the  equations  of  motion  will  be  the 
same  in  any  case.  The  horizontal  and  vertical  coordinates  of  a  point  on 
the  line  are  called  x  and  y,  respectively,  and  the  angle  between  the 
horizontal  and  the  tangent  to  the  line  is  designated  by  0.  Figure  2 
illustrates  the  corresponding  discrete  model  for  which  the  equations  will 
actually  be  derived.  The  line  is  divided  into  segments  in  such  a  way  that 
there  will  always  be  an  integral  number  of  them  between  any  points  where 
an  abrupt  change  occurs  in  some  parameter.  The  junctions  between  the 
segments  are  numbered  according  to  the  subscript  index  j,  which  runs 
from  0  at  the  anchor  to  s  at  the  surface. 

Before  we  can  properly  invoke  Newton's  law  of  motion,  it  is  necessary 
to  consider  the  inertial  properties  of  the  fluid  in  which  the  cable  is  immersed. 
We  shall  assume  that  the  kinetic  energy  of  the  surrounding  medium  is 
independent  of  the  component  of  velocity  parallel  to  the  line,  whereas  it 
varies  as  the  square  of  the  component  of  velocity  at  right  angles  to  the  line. 
Thus,  when  an  element  of  the  cable  is  accelerated  longitudinally,  no  hydro- 
dynamic  reaction  occurs,  but  when  the  cable  is  accelerated  laterally,  it 
behaves  as  though  it  possessed  additional  inertia.  The  component  of 
acceleration  normal  to  the  line  is 

anormal  =  -  x  sin  0  +  y  cos  0 
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The  accompanying  inertial  reaction  can  be  resolved  into  horizontal  and 
vertical  components.  Each  of  these  will  oe  proportional  to  the  corresponding 
component  of  aROrma^  namely, 

HORZ.  COMP:  -  anorma^  sin  0  =  x  sin^1  0  -  y  sin  0  cos  0 

VERT.  COMP:  anormal  cos  q  =  y  cos2  9  _  x  sin  0  cos  0 

Thus,  each  component  of  the  hydrodynamic  reaction  depends  on  both  components 
of  acceleration.  In  general,  the  reaction  force  is  not  parallel  to  the  acceleration 
vector  (except  when  the  tangential  component  is  zero),  so  that  it  is  necessary 
to  regard  the  inertial  parameters  of  the  system  as  tensors  rather  than  simple 
scalars. 

The  differential  equations  governing  the  motion  of  the  j station  on  the 
line  (see  Fig.  2)  can  be  written  in  matrix  notation  as  follows: 


where: 


8 


Each  lumped  mass,  mj,  has  been  expressed  as  the  average  mass  of  the  two 
segments  of  cable  which  lie  on  either  side  of  station  j .  Also,  one-half  the 
equivalent  transverse  mass,  ej,  of  the  entrained  fluid  associated  with  each  of 
these  segments  has  been  included  in  the  inertia  tensor.  Furthermore,  at 
those  stations  from  which  a  weight  is  suspended,  the  effective  horizontal  and 

v  Y 

vertical  masses,  m^  and  m-  ,  of  the  weight  are  to  be  added.  For  simplicity 
in  allowing  for  virtual  inertia,  we  have  assumed  that  any  such  weights  possess 
a  certain  degree  of  symmetry  and  remain  upright  as  the  line  moves  about. 

The  force  vector,  Fj,  on  the  right  side  of  eq.  (2. 1)  can  be  expressed 
as  the  sum  of  internal  forces  (the  tensions  acting  between  adjacent  mass 
elements)  plus  whatever  external  forces  are  present.  Thus,  in  expanded 
form  the  equations  of  motion  can  be  written 
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*i  *j  -  Vi  =  Tj  +  1/2  cos  ®j  + 1/2  -  Tj  -  1/2  cos  0j-  1/2  +  Xj 

(2.  2) 

-Kjiij+Jjjj  =  T j +  i/2  sin  ®j  +  1/2  -  Tj  -  1/2  sin  9j  -  1/2  +Yj 
where  Tj  +  -  tension  in  segment  of  line  between  stations  j  and  j  + 1 

Xj  =  horizontal  component  of  resultant  external  force  at  station  j 
Yj  =  vertical  component  of  resultant  external  force  at  station  j 
There  are  two  sources  of  external  force,  namely:  1)  gravity,  which  gives 
rise  to  the  weight  minus  the  buoyancy  and  acts  only  in  the  vertical,  and 
2)  fluid  resistance,  which  gives  rise  to  the  damping  forces.  Thus,  we  write 


where  Wj  =  img  -  1  p  g  (L  + 1  6-  1  +  1.1  d  1) 

J  2  J  2  2  J"2  ]~2 

Wj*  =  mj*  g  -  p  g  Vj* 

and  +  _1  =  drag  on  segment  of  line  between  stations  j  and  j  +  1 

2 

Xj*  =  horizontal  component  of  damping  force  on  weight  at  station  j 

Yj*  =  vertical  component  of  damping  force  on  weight  at  station  j 

Again,  in  order  to  get  the  best  approximation  to  the  continuous  case,  the  net 
effect  of  the  drag  at  station  j  has  been  expressed  as  one-half  of  each  component 
of  the  drag  on  the  segments  which  lie  on  either  side  of  this  station.  The 
buoyant  force  of  the  displaced  fluid  has  been  treated  likewise. 
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The  positive  normal  to  the  line  has  been  arbitrarily  taken  to  be  directed  upward 
when  0  equals  zero.  The  introduction  of  the  minus  sign  and  the  use  of  the 
absolute  value  of  one  of  the  velocity  factors  ensures  that  the  drag  will  always 
be  opposed  to  the  direction  of  qj  + 1/2  and  thus  act  as  a  dissipative  force  to 
remove  energy  from  the  system.  Since  the  velocities  of  the  two  endpoints  of 
each  segment  will,  in  general,  differ  slightly,  their  mean  value  (which  for  a 
straight  line  segment  is  exactly  equal  to  the  velocity  of  the  midpoint)  is  taken 
as  a  representative  value  in  the  definition  of  qj  +  \/2-  In  addition,  the 
definition  allows  for  the  presence  of  a  uniform  horizontal  current,  c,  to 
incorporate  the  ability  to  treat  towing  lines  as  well  as  mooring  lines  (or 
mooring  lines  subjected  to  ocean  currents). 
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In  addition  to  the  drag  on  the  line  itself,  there  will  also  be  resistance 
to  the  motion  of  any  concentrated  loads  which  may  be  suspended  from  the  line. 
These  additional  damping  forces  will  vary  with  the  velocity  but  will  not,  in 
general,  be  directed  exactly  opposite  to  the  motion  of  each  weight.  However, 
on  account  of  the  assumed  orientation  and  symmetry  of  any  such  weights,  the 
resistance  force  will  be  parallel  to  the  velocity  vector  whenever  the  relative 
motion  is  either  purely  horizontal  or  purely  vertical.  Accordingly,  the  two 
components  of  resistance  may  be  written 


Xj*  =  -  f  X  Uj  (5j  -  c) 

V  =-fjYVj 


where: 


fX  =  1  or.XcX 

fJ  2  PCj  T 
f.Y  _  1  q  p.Y  S.Y 

2  PC]  b] 

Uj  -  [  (ij  -  c)2  +  yj2  J  1/2 


(2.  5) 


Up  to  this  point  an  explicit  formula  has  been  given  for  the  evaluation 
of  every  term  in  the  equations  of  motion  (2.  2)  with  the  exception  of  the 
tensions.  To  determine  these  we  must  invoke  the  in  extensibility  condition 
which  was  assumed  at  the  outset.  This  takes  the  form  of  a  constraint  on  the 
motion  of  the  line.  It  requires  that  the  separation  between  adjacent  stations 
must  not  change  with  time.  Thus,  we  write 

(xj  -  x.  _l)2  +  (yj  -  y-j  _  i)2  =  1  j  _  1/2  =  const.  (2.6) 

This  holds  for  each  segment  of  the  line,  and  we  require  that  the  corresponding 
set  of  tensions,  Tj  _  take  on  values  such  that  the  resulting  solution  of 
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the  equations  of  motion  will  be  consistent  with  eq.  (2.  6).  Because  of  the 
implicit  nature  of  this  condition,  we  are  led  to  a  system  of  algebraic 
equations  for  the  determination  of  the  proper  tensions.  At  the  extremities 
of  the  line  (j  =  0  and  j  =  s)  xj  and  yj  must  be  obtained  from  the  boundary 
conditions,  namely: 


x0  =  x0(t) 

yo  =  yoW 

xg  =  xs(t) 

ys  =  ysW 

These  are  given  as  functions  of  time,  and  permit  the  introduction  of  any 
desired  types  of  driving  motions. 

Finally,  to  complete  the  formulation  of  the  problem  a  set  of  initial 
conditions  must  be  given  for  each  station  on  the  line.  Since  the  equations 
of  motion  are  of  the  second  order,  it  is  necessary  to  specify  both  the 
coordinates  and  the  velocities  at  t  =  0.  That  is, 

Xj  (0)  =  XjO  0  <  j  <  s 

yj  (0)  =  yf 

Xj  (0)  =  XjO 

yj  (0)  =  yj° 

where  the  superscript  index  "0"  is  used  to  designate  a  value  at  the  origin 
in  time. 


(2.  7) 


(2.8) 
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3.  SOLUTION  OF  EQUATIONS  BY  FINITE  DIFFERENCES 


( 


A.  General  Description  of  Computational  Procedure 

The  equations  governing  the  motion  of  a  cable,  as  derived  in  the 


last  section,  are  summarized  here.  The  basic  equations  of  motion,  equations 
(2.  2),  are  repeated  for  convenience, 


Ij  xj  -  Kj  yj  =  Tj  + 1  cos  0j  +  1  -  Tj  _  1  cos  0j  _  j.  +  Xj,  j  =  1,  2, . . .  S-l 

a  L  6  2 


(2.2) 


Kj  Xj  +  Jj  y.  =  Tj  +J.  sin  0j  +_1  -  Tj_ _1  sin  0j  _  1  +  Yj,  j  =  1,  2, ...  S-l 
2  2^^ 


where 


a)  S  is  the  number  of  junction  points 

b)  Ij,  Kj,  Jj  are  given  in  equation  (2.  1)  and  are  functions  of  the  physical 
properties  of  the  cable  and  of  position  only 

c)  Xj,  Yj  are  given  by  equations  (2.  3),  (2.  4), and  (2.  5)  and  are  functions 
of  the  physical  properties  of  the  cable  and  of  position  and  velocity. 

In  addition  the  motion  is  governed  by  the  condition  of  inextensibility  of  the 
cable,  equation  (2.  6), 

2 

(xj  -  Xj  _  j)2  +  (yj  -  yj  _  j  )2  =  Ij  _I  =  const,  j  =  1,  2, . . .  S.  (2.6) 

2i 

The  differentiated  (with  respect  to  time)  forms  of  this  relation 

(Xj  -  Xj  _  i)  (Xj  -  Xj  _  x)  +  (yj  -  yj  _  X)  (Yj  -  Yj  _  X)  =  o,  j  -  1,  2. . .  S  (3. 1) 

(Xj  -  Xj  _  j)  (Xj  -  Xj  _  i)  +  (yj  -  yj_j)  (yj  -  Yj  _  i)  +  -  Xj_  :)2  (3.  *) 

+  (yj  -  yj_i)2  =  °  j  =  l,  2. . .  s 

are  also  used  in  the  computation. 
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It  is  assumed  that  the  boundary  and  initial  conditions  are  known.  These  are 
given  in  equations  (2.  7)  and  (2.  8),  respectively.  The  system  of  equations 
summarized  above,  consisting  of  equations  (2.  2 ) ,  (2.  6),  (2.  7),  (2.  8),  (3. 1), 
(3.  2),  (3.  3).  (3.  4)  with  the  auxiliary  equations  (2. 1),  (2.  3),  (2.  4),  (2.  5) 
completely  describe  the  motion  of  the  cable. 


The  computational  procedure,  as  developed  in  detail  in  the  remainder 

of  this  section,  consists  of  an  algorithm  to  determine  the  values  XiH+1, 
11  1 
yjn  +  *,  x^n  +  2  ,  yin  +  ^  (at  time  t  =  tn+*  =  tn  +  A  t  and  tn  + 1  =  tn  + 

from  known  values  xj1,  yj1,  xj  2  yj  2  (at  time  t  =  tn  and  tn_  2). 

It  is  convenient  to  divide  this  algorithm  in  two  phases,  or  steps.  In  the  first 

step  tentative  (or  starting)  values  for  Xjn+*,  yjn+*,  Xj  2  ,  yj  2 

obtained.  In  the  second  step  improved  solutions  are  obtained. 


Step  1.  Using  equations  (2.  2)  and  (3.  2)  (3S-2  equations)  we  compute 

the  (3S-2)  unknown  variables  Xjn ,  y^n  (j  =  1,  2, . . .  S-l)  and  T^  l  (j  =  1,  2, . . .  S). 

2 

We  now  use  equations  (3.3)  and  (3.  4)  (4S-4  equations)  to  compute  the  (4S-4) 
variables  at  the  next  time  step  xjn+*,  yjn  ,  xjn+  2  ,  y j n  +  2  (j  =  T,  2. . .  S-l). 
These  are  considered  only  tentative  values  (denoted  in  subsequent  text  by 
use  of  the  tilde). 
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Step  2.  To  obtain  .the  improved  values  of  the  tensions  Tf_  1  (j  =  1,  2, . . .  S)  I 

1  1  J“2 

and  the  quantities  xjn- ,  yjn  ,  Xjn+*,  yjn+*,  Xjn+2,  yjn+2^  (a  total  of  (7S- 6) 
quantities)  we  use  the  system  of  equations  (2.  2),  (2.  6)  and  (3.  3),  (3.  4), 
consisting  of  (7S-6)  equations.  Hovever,  since  equations  (2.  6)  are  not  linear 
but  quadratic  in  the  unknowns  Xj  yj  an  explicit  solution  is  impractical  to  obtain. 

For  this  reason  a  computation  algorithm  based  on  the  Newton- Raphson  method 
of  successive  approximations  is  developed.  A  detailed  discussion  of  the  compu¬ 
tation  procedure  used  in  this  problem  is  given  in  the  sections  which  follow. 

B.  Determination  of  Tentative  Values  of  Tensions 

In  summary,  the  method  of  solution  at  each  time  step  involves  in  the 
first  phase,  1)  the  determination  of  a  tentative  (but  consistent)  set  of  tensions 
Tj  - 1/2  for  segments,  and  2)  numerical  integration  of  the  equations  of 
motion  to  predict  xj  and  yj  one  step  ahead;  in  the  second  phase  3)  evaluation 
of  the  discrepancies  in  the  constraint  equations  from  which  a  set  of  first 
order  corrections  to  the  tensions  can  be  obtained,  and  4)  integration  of  the 
equations  a  second  time  to  obtain  corrected  values  of  the  coordinates. 

The  system  of  equations  (2.  2)  may  be  regarded  as  a  set  of  (2S-2) 
linear  equations  in  the  variables  xj,  yj  (accelerations)  and  may  be  solved 
directly  for  these  variables.  If  we  designate 
Lj  =  (A  t)2  Ij  /  (Ij  Jj  -  Kj2) 

Mj  =  (At)2  Jj  /(Ij  Jj  -  Kj2) 

Nj  -  (At)2  Kj  /  (Ij  Jj  -  Kj2); 
then  the  equations  of  motion  (2.  2)  can  be  reduced  to: 
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*j  =  [Rj  TJ  + 1/2  -  PjTj-t/2  +  Uj]  /  (A t)2 
h  ■  [Sj  Tj+1/2  -  QjTj.t/j  +  Vj]/(4t)* 
where:  Pj  =  Mj  cos  9j-i/2  +  NjSin8j_i/2 

Qj  =  Nj  cos  0j- 1/2  +  Lj  sin  9j_  1/2 

Rj  =  Mj  cos  0j+i/2  +  Nj  sin  0j+1/2 

Sj  =  NjCosGj+1/2  +  LjSin0j+1/2 
Uj  -  Mj  Xj  -  Nj  Yj 

vi  =  Nj  Xj  f  Lj  Yj 

We  observe  that  equation  (3.  2)  involves  positions,  velocities,  and 
accelerations.  As  is  often  the  case  with  finite  difference  procedures,  it 
proves  to  be  convenient  to  compute  positions  and  accelerations  at  the  mesh 
points  while  velocities  are  found  at  the  mid-points  in  time.  For  this  reason 
we  shall  use  a  modified  form  obtained  by  evaluating  this  equation-  at  t  =  tn 
and  t  =  tn_*,  and  then  adding  the  two  results  together,  namely, 


,  n  n  \/-n  ~n  x  ,  n  n  w..  n  ..n  n-1  n-1  n-1  ..n-l, 

(Xj  -  Xj_i) (Xj  -  Xj_x)  +  (yj  -yj_i)(yj  -yj-i)  +  (Kj  -  xj-i  )  -*j-i) 

+  (yj1  1  -  yj1- 1)  (yj1-1  -  yj1-'/ )  +  2 (At)-2  [ (xj1  -  xj1-1)  -  (xj1.!  -  xj1;/)] 

+  2 (At)" 2  [(yj1  -  yjn_1)  -  (yj1-!  -  yj1-/)]2  =  0;  (3.6) 

in  which  we  have  used  the  approximations, 

(ij'-Xj.i)2  +  (Xj"  1  -  x"./ )2  =  2(x“"  2  -  x".!2)2 


(y"  -  v?1  ,)2  +  ft"'1  -  v"',1)2 
yj“l'  'yj 


„  r,  n  n-1,  ,  ,  n  n-1  .  1 2 

2  [(xj  -  Xj  )  /At  -  (xj_1  -  Xj_j  )  /At] 

n-1  n-1  9 

2(yj  ^  -  y j _ ! zr 


■  2  [(yj*  -  yf  't/it  -  (y"_ !-  y".",1  )/At]2  . 


n  n-1 


Note  that  equations  (3.  6)  are  linear  in  the  accelerations.  Likewise  equations 
(3.  5)  are  linear  in  the  tensions.  Consequently,  when  these  expressions  are 
substituted  for  the  acceleration  components  in  the  constraint  equations  (3.  6), 
we  obtain  a  set  of  conditions  which  are  linear  in  the  tensions,  namely, 


E: 


,n  ,  -11  Tn  4.  TP11"1  1 

■j-l/2  Tj-3/2  '  Fj-1/2  Tj-l/2  +  Gj-l/2Tj+1/2  Ej-l/2  Tj-3/2 


'  ^-1/2  Tjn-"l/2  +  GH/2  Tj1+"l/2+  H?- 1/2  +  Hjn-’l/2  +  2^"  -  xj1'1)  -  (x"-l  -  xj1;/)]2 


„  r.  n  n-L  /  n  n-Ll  2  _ 

+  2  Ltyj  -  yj  )  -  (yj-i-yj-i)J  =  o 

where:  e"_1/2  =  (xj1  -  xj1.,)  p"_,  +  (yj1  -  yjL )  Qj-i 

_,n  .  n  n  .  ,„n  n  .  ,  n  n  ,  n  n  , 

Fi-i/2  =  <Xj  -Xj-iXPj  +hj-i)  +  (yj  -  yj-j)  (Qj  +  Sj.,) 


(3.  7) 


n 

Oj-l/2 


■  (xf  -  x?  1)  R?  ♦  (yj1  -  vU  $ 


<-1/2  =  (x"  -  xj1.!)  (of  -  Uf.j)  f  (y"  -  yj1.,)  (<  -  <_,) 

Now  assume  that  the  solution  is  correct  up  to  t  =  tn.  Then  all  quantities 
in  (3.  7)  can  be  evaluated  at  once  except  for  Tj1.  3/2  T j1, 1/2  and  ^j+z/2  • 

The  tentative  values  of  the  tensions  -  signified  by  the  tildes  —  are 
determined  by  the  following  system  of  equations: 


f  n 
-  F, 


n 


0.  5 


G0.  5 


E1!.  5  "  *1.  5  Gl.  5 


En  =  -  F 


n  ~n 


2.  5  '  r2.  5  u2,  5 


/Tn  \ 
l0.5 


Tl.  5 
n 


t2.5 


/  n 
-*0.5 


-  T 


n 

1.5 


Eg_i_  g  -  F^_j  g  Gs_i  g 


Es-0.  5  "  Fs-0.  5 


n 


- 


*  n 
*  2.  5 


.n 


rp11 

Ts-1.  5 


,n 


VTs-0.5j 


>  0.8) 


s-1.5 

s-0.5 
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where  S V 


n  n-1  J1"1  Tn-1  ,  rn-l  Tn-1  „n-l 

j-1/2  =  Ej - 1/ 2  Tj - 3/ 2  Fj-l/2  Tj -1/2  Gj-l/2Tj+l/2  j-l/2 

„n  .n  n  n-1.  .  n  n-1  1 2  r  n  n-1  n  n-li2 

nj-i/2 +  2n  -  xj  !  -  (xh ■  xj-i >]  +  2  h  -  yj  )  -  frj-i-  yj.pJ 


In  general,  we  can  write:  (fori  j  fC  s) 
Ef-l/2  Tj-3/2  -  rf-1/2  tf-1/2  +  g“-1/2  Tj" 


+1/2  +  ^  j-1/2  “  0  <3-9) 


n  n 

with  the  conditions:  Eg  5  =  Gs_j/2  =  0  for  all  n 
Also  Pg,  Qq,  Rq,  So  andPg,  Qg,  Rg,  Sg  =  0  for  all  n. 
and  Ug  =  (At)^  Xq  and  =  (At)2  Xg  for  all  n. 

Vq  =  (At)2  yg  and  Vg  =  (At)2  yg  for  all  n. 

The  matrix  of  coefficients  of  the  system  of  equations  for  *s  a  triple 

diagonal  one,  and  it  can  be  easily  reduced  to  a  single  linear  equation  by 
elimination.  Thus,  we  solve  equation  (3.  9)  for  Tj+j/2  . 


_  J-*/2  v71  / 

‘j+l/2  "  „n  ‘j-1/2 

GJ - 1/ 2 


P11  T  11 

Ej ■  1/2  *n  *  j-1/2 

rn  '  -  j-3/2  pn 

Gj  - 1/2  j-1/2 


(3. 10) 


Now  we  express  each  tension  as  a  linear  function  of  Tg  5  (the  tension  in  the 


first  link)  as  follows: 


1+1/2 


=  Of. 


Tn  c 


j+l/2  A0. 5  +  3+I/ 2 


~n  n 

T  a  c  +  Bi 


.j-1/2  -  “j-1/2  V/T0. 5  +  “j-1/2 

~n  n  ~n  n 

Tj -3/ 2  "  “j-3/2  T0. 5  +  fij-3/2 


(3.11) 


and  we  arrive  at  the  following  recursion  formulas  for  «j1+i/2  anc*  ^j+l/2>  namely, 
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n  /Tn  n  n  n  .  ,  n 

“j+t/2  =.  (Fj-l/2  “j- 1/2  -  Ej-1/2  aj-3/2)/Gj-l/2 


(3. 12) 


<W  =  (f"-1/2  ^.l/2-EU/2^/2-^-l/2)/GU/2 

with  the  conditions:  ag.  5  =  1,  a_gt  5  =  0  for  all  n. 

n  n 

fig.  5  =  0  fi_o.  5  =  0  for  all  n. 

Starting  with  j  =  1,  we  evaluate^  OjVi/2  ana  tf+l/2  recursively  up  to  j  =  s  -  1. 
~n 

We  then  find  To.  5  from  the  last  equation  of  the  system,  (j  =  s),  using  the  same 

substitutions  as  before,  that  is, 

~n  n  ~n  „n 

Ts-3/2  “  aj-3/2  T0. 5  +fij-3/2 


~n 

T 


n 


;n 


n 


s-1/2  =  aj-l/2  T0. 5  +fij-l/2 


The  final  result  is 

~n  (FS-l/2  ^s- 1/2  '  Es-1/2  Bs-3/2  "  * 

T0. 5  ~  - 


s-1/2) 


n  n  „n  n  . 

(fs-1/2  <*s-l/2  "  Es-l/2«s-3/2  ) 


(3.  IS) 
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C.  Method  of  Determining  Improved  Tensions 


In  order  to  solve  eqs.  (3.  5)  numerically,  we  replace  xj1  and  yj1  by 
their  simplest  central  difference  approximations  (equation  3.  4),  namely, 


Xj  =  -  2XJ1  +  X-1  A)  /  (At)2 

-n  .  n+1  ,  n  1-1.  .  <*■ 

^  =  “  2yj  +  yj  )/  (At) 

Now  we  solve  for  Xj1**  and  yj1** ,  considering  these  as  tentative  values  subject 
to  a  slight  modification  in  order  to  satisfy  a  system  of  constraints.  Thus  we 
write 


n+1 


n-1. 


_n+l  n  n-1 

Xt  =  2x4  -  X4 


1 


1 

n 


3 


n  ~n  n  ~n  n 

Pj  Tj-l/2  +  Rj  Tj  +1/2  +  Uj 


(3. 15) 


n  ~n 


n  ~n 


n 


.n+1  n  n-1  ..  ..  „  __ 

yj  =  2yj  -  yj  -  Qj  TjA/2  +  Sj  TJ+1/2  +  V, 

n  n  n  n  n  n 

The  quantities  Pi ,  Q; ,  R. ,  SC ,  IL  and  V.  are  the  same  as  were  used  to 

J  J  J  J  J  J 

set  up  the  coefficient  matrix  for  the  tensions,  and  the  values  for  and 

fj+1/2  are  obtained  from  eqs.  (3. 11). 

Next,  we  determine  the  set  of  corrections  6 T^_  1/2  to  be  applied  to  the 

~n  n+l  n+1 

tensions  Tj _  1/2  in  order  that  the  values  of  xj  and  yj  should  also  satisfy 


the  inextensibility  condition  (2.  6).  For  this  purpose  we  define  the  function 
n+1  1  r  n+1  n+1  2  ,  a+l  n+1  2  2 

nj.l/2  *  2  [(xi  -XH>  +  ■  yM>  -  *1-1/2. 


(3. 16) 


which  measures  the  discrepancy  in  the  distance  between  the  extrapolated 

positions  of  pairs  of  adjacent  stations.  We  observe  from  eqs.  (3. 15)  — 

n+1  n+1 

with  the  tildes  suppressed  —  that  Xj  and  yj  are  functions  of  the  tensions, 


that  is, 
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( 


n+1  n+1  /  n 

xj  -  xj  tTH/2 . 

n+1  n+1  f  n 

yj  -  yj  {Tj-i/2>  Tj+i/2/ 

n+1 

Consequently,  -^j-l/2  may  a^s0  t>e  expressed  as  a  function  of  the  tensions. 
This  enables  us  to  write  the  system  of  constraints  to  which  the  tensions  are 


TH/2) 

n 


subject  as  follows: 


17. 


n+1 


j- 1/2 
n+1 


,  n+1 
“1-1/2 


TH/2>  T"l/2-  T"+-l/2 }  *  0  ^  P-17> 


n+i 

since  17  .  j /£  vanishes  when  the  inextensibility  condition  is  obeyed. 


Now  let  T, 


H 

V  3 

CO 

to 

It 

~n 

Tj - 3/2  + 

,  rr.n 

6  Tj-3/2 

n 

» 

-n 

.  n 

Tj-l/2  = 

Tj-l/2  + 

6Tj-l/2 

n 

~n 

k  mn 

Ti+l/2  ’ 

Tj+l/2  + 

5  Tj  +1/2 

(3. 18) 


n+1  r~  n  ~n  ~n  ) 

:xpand  17  j_2/2  ^  a  Taylor  series  about  the  point  | Tj_3/2 »  Tj - 1/2»  Tj+l/2 j- 

,  we  obtain  ~n+*  ~n+* 

+1  .n+1  8^,i-l/2  n  9%l/2  n 

-1/2  =  ^ j - 1/2  3Tjn_3/2  5Tj"3/2+  6TH/2 


j-1/2 


...  n+1 


(3. 19) 


+ 


9T? 


h,! 


-  n+1  n+1 

^ j-1/2  s  ^ 


in 

j+1/2 

higher  order  terms. 

n+*  f  ~n  ~n  ~n  -t 

j-1/2  |Tj-3/2i  Tj - 1/2  >  Tj+l/2  | 

1  i\~n+l  -n+lv2  rn+l  ~n+lN2  ,  2  ] 

=  i  l(xj  -xi-i>  +  <yj  -yj-P  -  ij-vzJ 
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Provided  that  the  tentative  values  Tj!_j/2  are  sufficiently  close  to  the  correct 

values  Tjij/2  j  we  may  neglect  the  higher  order  terms  in  the  expansion  (3. 19), 

and  thereby  obtain  a  system  of  s  linear  equations  for  the  differential  corrections 

5  T^-i/2  •  These  equations  have  the  same  form  as  the  previous  system  (3.  8) 
~n 

for  determining  namely, 


-n+l 
“F0.  5 

-n+l 
G0.  5 

5T0. 5  1 

i 

-n+l  ! 
-%  5 

~n+l 
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1 
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1 
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the  general  expression  being:  (for  1  <  j  <  s) 


-~n+l  .  n 

Ej-l/2  6Tj-3/2 


~n+l  n  ~n+l  n  -n+l 

Fj-l/2  5Tj-i/2  +  Gj-1/2  6Tj+i/2  +  ftj-1/2  =  0 


(3.  21) 


where : 


~n+l 


3Tj-3/2 


#-n+l  -n+L 

(Xj  -  Xj.j) 


_n+l 

_ 

aTH/2 


.n+l 

exi-l 

3T?-3/2 


,~n+l  ~n+L 

+  (yj  -yH) 


3y 


n+l 

j 


3f-3/2 


_n+l 

ayj-i 

3T"-3/2 


0 
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~n+l  ..n+1  _n+l  .  n  n  n+1  -n+1.  ,  n  n 

Fj-l/2  =  (xj  “xj-l^Pj  +Rj-l)  +  (yj  -  yj-1 )  (Qj  +Sj-l) 


(3.  22) 


-n+1  .n+1  .n+1  n  n+1  _n+l  n 

Gj-l/2  =  (Xj  -  Xj_i)Rj  +  C/j  -Yj-l  >Sj 

-n+1  _n+l 

with  the  conditions:  Eg,  5  =  Gs-l/2  =  0  for  all  n,  and  the  quantities 
pj1,  Q^,  and  sj1  being  the  same  as  before. 

The  system  (3.  20)  can  be  solved  in  a  manner  completely  analogous  to 
the  solution  of  the  system  (3.8).  Thus,  we  write 

STH/2  =  'j-1/2  6T0.5  +  X"-l/2  P+3> 

and  obtain  the  following  recursion  formulas: 


24 ' 


(3.  24) 

n  ~n+l  n  ~n+l  n  =.n+l  ~n+l 

Aj+l/2  “  (Fj-l/2  xj-l/2  ~  Ej - 1/ 2  xj-3/2  "  ^j-l/2)/Gj-l/2 

with  the  conditions:  kq.  5  =  1,  *-0.  5  =  0  for  all  n- 

*0.  5  =0.  »-0. 5=0  . 

Finally,  the  last  equation  of  the  system  (when  j  =  s)  enables  us  to  solve  for  6T(j,  5. 
The  result  is 
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For  solution  on  an  automatic  computer  it  is  more  convenient  to  express 
eqs.  (3.  27)  in  terms  of  corrections  to  be  added  to  the  tentative  values  of 

the  coordinates.  That  is, 

6  xn+l  =  .  pH  5  -I'll  +  Rn  gTn 

j  j  j-1/2  j  j+1/2 


5  yn+1  =  -  Qn  5  Tn  . 
y3  1  j+1/2 


+  S11  6t 
j 


n 

j+1/2 


Then  the  corrected  coordinates  are  given  by: 


xn+l  = 

jn+l  +  6xn+l 

i 

J 

1 

(3.  29) 

yn+1  - 

yn+1  ♦  «yn+1 

1 

J 

J 
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These  values  are  now  accepted  as  final.  Moreover,  as  soon  as  the  values 

n+j 

of  ,  (to  be  used  with  eqs.  (3. 8)  for  the  next  time  step)  have  been 

j  - 1/2 

computed  and  stored,  the  cycle  of  computations  is  finished  and  there  is 


no  further  need  to  retain  the  values  of  P*1, 

J 


Q?,  Rn,  S11  and  T 
1  J  3 


n 

j-l/2* 


9 

r 


E.  Special  Form  of  Equations  for  Computing  First  Time  Step. 


We  assume  that  the  velocity  components  are  zero  at  each  station, 
and  obtain  the  initial  coordinates  from  the  equations  for  static  equilibrium 
of  the  line.  Since  x?  and  y?  =  0,  eq.  (3.  2)  reduced  to 

(xj  •  xj-l>  (xi  -  xj- 1;'  -  <yj  -  yj-l>  (yj  -  yj-l> =  0 

and,  on  substituting  the  expressions  (3.  5),  we  find  that  the  tensions  are 
subject  to  the  constraint 


(3.30) 


pO  rpO  ■pO  itiO  ,  rp  np®  4.  tj® 

J-l/2  j-3/2  j-l/2  Tj- 1/ 2  3 -1/2  j+1/2'  j-l/2 


=  0 


(3.31) 


Comparing  this  with  eq.  (3.  9),  we  see  that 


i  0  =  H" 

j-l/2  j-l/2 


To 


The  system  of  equations  (3.  8)  is  then  solved  in  the  usual  way  to  get  the 
proper  initial  tensions  Tj 


(3.32)  f  1) 


■J 

To  obtain  tentative  values  for  the  coordinates  at  t  =  tA ,  we  make  use  of 


their  Taylor  series  expansions  about  the  point  t  =  t°,  namely: 


x?  =  X?  +  (At)  X?  +  i(At)2x°  +  ... 

y \  -  y°  +  (At) y?  +  I(At)2y°  +  ... 

J  J  J  9  J 


Taking  x? and  y°  =  0,  and  substituting  eqs.  (3.  5)  for  x°  and  y?,  we  find 
J  j  J  J 


(3.3*) 
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<*> 

u 


xl  = 


x°  +  * 
J  2 

y°  +  — 
2 


PJ°Tj°-1/2  +  Rj)Tj>+l/2  +  U°. 


Qj°  Tj  - 1/ 2  +  S°]TU/2  +VJCJ 


(3.  S*) 


The  corrections  to  the  tensions  are  then  determined  by  the  system  of 
equations  (3.  20)  in  the  usual  manner.  Finally,  the  corrections  to  the 
coordinates  are  computed  as  follows: 


5x!  = 
] 

«yj  = 


1 

ft 

L 

1 

2 


P?  6Tj0-l/2 


+  R°  6T°  1  /9 
1  3+1/2 


/  6TH/2  +  Sj  %°+l/2 


(3.35) 


and  the  corrected  coordinates  are  given  by: 


x.l  = 


Yi 


1  _ 


x}  +  6x} 


y}  +6y1 

J  i 


(3.  36) 
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ANALYSIS  OF  NUMERICAL  STABILITY 
In  order  to  obtain  a  valid  solution  of  the  system  of  partial  differential 
equations  (2.  2),  (2.  6),  governing  the  generalized  motion  of  a  cable  it  is 
necessary  to  insure  the  staoility  (in  the  sense  discussed  in  References  1,  2,  3)* 
of  the  equivalent  finite  difference  system  (2.  2),  (2.  6),  (3.  3),  (3.  4).  In  this 
section  we  will  derive  the  criteria  for  stability  of  this  system  of  equations. 

We  will  also  show  that  whereas  the  system  of  finite  difference  equations 
(2.  2),  (2.  6),  (3. 3),  (3. 4)  is  stable  for  sufficiently  small  time  intervals  At, 
the  system  (2.  2),  (3.  2),  (3.  3),  (3.  4)  is  always  unstable.  This  characteristic 
of  the  latter  system  has  led  to  the  abandonment  of  this  simpler  set  of  equations 
in  favor  of  the  more  difficult  nonlinear  system  (2.  2),  (2.  6),  (3.  3),  (3.  4). 

In  order  to  determine  the  stability  of  a  system  of  finite  difference 
equations  we  study  the  growth  of  a  small  disturbance  or  perturbation.  The 
conditions  for  stability  are  said  to  be  satisfied  if  the  amplitude  of  a  small 
disturbance,  introduced  at  any  time,  t,  in  any  of  the  dependent  variables, 
does  not  increase  exponentially  with  successive  time  steps.  This  condition 
may  be  stated  as  follows: 

If  5F(s,  t)  and  6F(s,  t  +  At)  are  values  of  a  variation  (or  perturbation) 
in  any  of  the  dependent  variables  x,  y,  T  in  the  system,  then  it  is  said 
to  be  stable  provided  |s  F(s,  t  +  At)/6F(s,  t)  <  1.  We  introduce 
perturbations  6x,  6y,  <5T  in  the  independent  variables  x,  y,  T,  respectively. 

For  the  sake  of  the  stability  investigation  we  further  assume  that  ej  is 
negligible  compared  to  mj.  Substituting  in  equations  (2.  2),  (2.  6),  (3.3),  and  (3.4) 
we  obtain  the  variational  system  of  equations 
*  References  are  listed  on  page  41. 
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)  ffij6*j  =  Tj+V2  6C0S  Vl/2  '  Tj-l/2  6C0S  Vl/2  +  C0S  8j+1/2  «Tj+1/2 

-  cos0j-l/2  6Tj-l/2  '  |[Dj+l/2«sin8i+l/2  +  Dj.i/2  6sin0j.1/2 
+  sin  9j  +1/2  6Dj+1/2  +  sin  l/2S Dj- 1/2] 

“j6yj  =  Tj+l/'26sin  0j+l/2  "  Tj.j/gesin  0j_1/2  +  sin  0j+l/26Tj+1/2 

(4. 1) 

-  sin  0j_1/2  5Tj_i/2  +  \  [Dj+l/26cos  ej+l/2  +  Dj-l/26cos  0j_1/2 
+  cos  0j+i/26Dj+1/2  +  cos  0j_i/25E>j_i/2] 

cos0j+l/26cos0j+l/2  +  sin9j+l/26sin9j+l/2  =  0 
where, 

6Dj+l/2  =  "  2fj+ 1/2  | ^j+1/2 1  6(*j+l/2 

=  mj  +  m j* 

Wj+1/2  =  -  f  [*Vl  "  c)  +  (V  c)]  6sin  0j+l/2  +|  &J+1  +  VScos  0j+1/2 
"  |  sin  0j+l/2  (6xj+i  +  5  xj)  +  I  cos  0j+1/2  (6y.+1  +  6y.); 

and  where, 

5cos  0j+1/2  =  (6xj+1  -  6xj)/lj+1/2,  5  sin  0j+1/2  =  (6yJ+1-  6yj)Aj+i/2; 

5x?-1/2  =  (6x?-  6x?_1)/At,  5yn_1//2  =  (5yn-  6yn_1)/At; 
j  J  J  J  J  J 

6x?  =  25xn  +  5xn_1)/(At)2,  5yn=(6yn+1-  2  5yn  +  5yn_1)/(At)2. 

J  J  J  J  J  J  j  j 
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We  will  assume  in  this  analysis  that  within  a  small  region  in  the  (s,  t) 


plane  the  coefficients  (T?,  cos  0j\  D?,  etc. )  of  the  variational  functions  vary 

3  3  3 

only  slightly  and  hence  may  be  treated  as  constants.  We  will  denote  these 

simply  by  T,  cos  0,  D,  etc. ,  omitting  the  subscripts.  A  solution  of  the 

system  of  equations  (4. 1)  can  then  be  obtained  in  the  form, 

6xn  =  aeifi3+anAt 
3 

6  n  =  beiBj+®iAt 

i 

BTj1  =  ceifii+tmAt 

where,  a,  b,  c  are  real  constants  and  a  complex.  Substituting  in  equation  (4. 1) 
we  obtain  a  system  of  linear  homogeneous  equations  for  the  quantities  a,  b  and  c 
which  has  a  non-trival  solution  provided  the  determinant  of  the  coefficients  is 
identically  zero.  After  some  algebraic  simplifications  the  determinant  of  the 


coefficients  may  be  written  in  the  form 


F  -  A  sin  9  D'  +  B  sin  0  cos  0 

-D'  +  A  cos  9  F  -  B  cos  0  sin  0  =0  (4.2) 

cos  0  sin  0  0 


where 

A  =  f  !q|  |^2 i y  sin  fi  -  (1  sin  9/At)(l  +  cos  fi)(l  -  A"*) 

B  =  f  [q|  1 2  i  (x  -  c)  sin  fi  -  (1  cos  0/At)  (1  +  cos  6)  (1  -  A”  *)  j 
D'  =  i  D  sin  fi 

F  =  mtS/(At)2  +  4T  sin^  ^ 

Lt 
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and  where 

X  =  eaAt,  |  =  (A.-  2  +  X-1). 

Multiplying  the  elements  of  the  determinant  and  simplifying  we  obtain 

A  sin  9  +  B  cos  0  -  F  =  0  (4.  3) 

But,  A  sin  G  +  B  cos  0  = 

f  |q,  j(2  i  sin  B)p  -  (i/At)(l  +  cos  6)(1  -  A'*) 
where  p  =  (x  -  c)  cos  0  +  y  sin  0 

i.  e. ,  the  tangential  component  of  the  velocity  of  the  cable  (relative  to  the 
medium). 

Substituting  in  equation  (4.  3)  we  finally  obtain 

ml  A2  +  /f  |q  At  [4(1  +  cos  B)  -  p  At  (2i  sin  B)j 

..  (4. 4) 

+  4T  (At)2  sin2  J*  -  2mi  j  X  +  Jmi  -  f  |qj  4  At  (1  +  cos  B)J  =  0. 

2 


Now,  comparing  the  first  and  second  terms  of  the  coefficient  of  X  we  find 
that  the  second  term  is  negligibly  small  provided  2p  At  «  i,  i.  e. ,  the 
tangential  distance  traversed  by  the  cable  in  one  time  step  is  very  small 
compared  with  the  length  of  the  cable  segment.  Since  this  is  usually  the 
case  and,  at  any  rate,  can  always  be  satisfied  by  taking  the  time  step 
sufficiently  small, we  will  omit  this  term  from  our  suosequent  analysis. 

For  the  case  of  negligible  drag,  i.  e. ,  f  =  0,  approximately,  we  obtain 
from  equation  (4.  4) 


X2  +  i4T(sin2  |)(At)2/ml  -  2  ]  X  +  1  =  0. 

Li 

In  order  for  the  solution  to  be  stable,  the  conditions  |Ajj  <  1, 
must  both  be  satisfied.  But  if  Xj  is  a  solution  of  (4.  5)  then  X2  = 


W  ^  1 

1  is  also 
X1 


(4.  5) 
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a  solution.  It  follows  that  the  conditions  for  stability  can  only  be  satisfied  if 

=  —  =  I Xi|  =  1.  Now,  let  Aj  =  cos  y  +  i  sin  y,  X2  =  cos  y  -  i  sin  y  =  t-; 


i.  e. ,  Aj  +  A£  =  2  cos  y  <2.  Again  from  equation  (4.  5) 


4  T  (sin2*h(At)2 


Aj  +  A2  -  2  - 


We  thus  obtain  the  inequality 


2  ~  4T  (sin2|)(At)2/mlr  <  2 
z 


T(sin^  -)(At)2/mi  <  1. 
2 


(4.  6) 


This  requirement  is  tantamount  to  the  condition, 


At  < 


velocity  of  transverse  wave  ' 


In  the  more  general  case,  allowing  for  finite  drag,  equation  (4.  4)  may 
be  reduced  to  (after  neglecting  the  second  term  of  the  coefficient  of  A), 
mix2  +  [f|q|lAt(l  +  cos  B)  +  4 T(At)2  sin2  ^  _  2ml]  A 

+  [ml  -  f  jq| i  At  (1  +  cos  B)]  =  0 

This  equation  is  more  difficult  to  analyze.  However,  it  is  possible  to  show 
that  both  Aj  <  1  and  A2 1  <  1  provided  the  slightly  more  stringent 


conditions 


At  <  JE*  ,  and 


At  < 


(4.8) 


2f|q| 


are  satisfied. 
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We  will  now  show  that  the  replacement  of  equation  (2.  6)  by  its 
differentiated  form  (3.  2)  results  in  an  unstable  system;  and  that  furthermore, 
the  use  of  any  time  interval  At  no  matter  how  small  does  not  change  the 
unstable  character  of  the  equations.  It  will  suffice  to  show  that  this  condition 
exists  in  the  case  when  the  drag  is  negligible,  i.  e. ,  f  =  0.  The  variational 
equation  corresponding  to  equation  (3.  2)  is, 

(xj  -  Xj_1)(6xj  -  6xj_j)  +  (xj  -  Xj_1)(6xj  -  Qxj.j)  +  2(xj  -  Xj.jHSxj  -  Sxj.j) 

+  (yj  -  yj_i)(5yj  -  tyj-i)  +  (yj  -  yj-iKfyj  -  Sy^)  +  2(yj  -  yj_i)(6yj  -  5yj_i)=  o. 


<sv 

'& 


Substituting  appropriate  values  for  5x,  5y  and  neglecting  terms  containing  f 
the  determinant  equation  (4.  2)  is  replaced  by 


0 


cos  0 


(xj  -  Xj„!)|  +  (xj  -  xj.jHAtr 
+  2(xj  -  Xj_1)(l  -  A-1)  At 


F  sin  9 

(yj  -  yj_i)^  +  (yj  -  yj-i)(^t)2  o 

+  2(yj  -  yj_i)(l  -  X'1)  At 


Multiplying  the  elements  of  the  determinant  we  obtain 
F  cos  0  [(xj  -  Xj_|)|  +  (Xj  -  Xj_  j)(At)2  +  2(xj  -  xj_i)(l  -  X_1)Atj 
+  F  sin  0  [(yj  -  yj_j)t  +  (yj  -  yj_i)(At)2  +  2(yj  -  yj_i)(l  -  A”1) At' 


=  0 


Equating 


cos  0  = 


a  _  yr  yj-i . 


,  sin  0  =— l - 1— L  ; 

h- 1/2 


and  using  the  relation  (first  time  derivative  of  equation  (2.  6)), 
(Xj  -  xj_ j)(xj  -  Xj.j)  +  (yj  -  yj.jKyj  -  y^)  =  0 


=  0 


(4.  2) 


(4.9) 
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as  well  as  equation  (3.  2)  we  obtain  in  place  of  equation  (4.  9) 


F  [l2S  -  [(xj  -  xj_j)2  +  (yj  -  yj_i)2]  (At)2  ]  =  0. 


Thus,  in  order  to  satisfy  the  stability  conditions  the  following  two  equations 


must  be  satisfied 


F  =  0 


(4. 10) 


12$  -  [(xj  -  xH)2  +  -  yj_!)2](At)2  =  0. 


(4. 11) 


It  can  be  shown  that  equation  (4. 10)  is  equivalent  to  the  criterion  (4.  6)  and  is 
satisfied  provided 


However,  equation  (4. 11)  can  never  be  satisfied  for  any  finite  At,  since  it 
requires  that 

£  =  [(Xj  -  Xj-l)2  +  (yj  -  yj_1)2](At)2/l2, 

a  positive  quantity.  This  conclusion  follows  as  a  result  of  the  definition 
£=X-2  +  A"*.  If  Ai  is  a  root  of  equation  (4. 11),  then  X2  =  is  also  a 
root  of  this  equation.  As  before,  it  follows  that  for  stability  Xj  |  <  1  and 


|  X2 1  =  y-  <  1.  Hence  jXj[  =  j  X2I  =  1.  Let  Xj  =  cos  7  +  i  sin  y, 

y2  =  cos  7  -  sin  7  =  ;  then  %  =  2  (cos  7  -  1),  or  -  4  <  %  <  0.  Thus, 

c  A1 

to  satisfy  the  stability  requirement  %  must  lie  between  0  and  -  4,  and 
consequently  is  always  negative  or  zero. 

In  Figure  3  the  vertical  velocity  of  the  midpoint  of  a  mooring  line  is 
plotted  as  a  function  of  time,  both  as  obtained  by  the  use  of  the  stable  (valid) 
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system  of  equations  (2.  2),  (2. 6),  (3. 3),  (3. 4)  and  as  obtained  on  the  basis 
of  the  unstable  (invalid)  system  (2.  2),  (3.  2),  (3.  3),  (3.  4).  It  will  be 
noticed  that  at  approximately  18  seconds  the  unstable  solution  rapidly 
increases  beyond  any  reasonable  limit. 


5.  RESULTS  AND  CONCLUSIONS 


A  number  of  solutions  were  carried  out  for  varying  wave  heights  and 
periods.  Several  typical  solutions  are  reproduced  here  for  the  information 
of  the  reader.  In  Figures  4,  5,  and  6,  plots  are  given  of  the  maximum 
tension  attained  along  the  cable  as  a  function  of  time  for  'ave  heights  of 
6  feet  and  periods  of  12.  5  seconds,  7.  5  seconds,  and  5  seconds,  respectively. 
The  periods  of  the  variation  in  maximum  tension  correspond  to  the  periods 
of  the  forced  vibration,  as  expected.  The  maximum  tension,  however, 
increases  in  amplitude  from  32,  250  lb  in  the  case  of  the  12.  5  sec  period 
waves  to  38,  500  lb  for  7.  5  sec  period  waves  to  49,  500  lb  when  the  period 
is  5  seconds,  hi  Figure  7  the  maximum  tension  attained  for  wave  heights 
of  9  ft  and  a  period  of  7.  5  seconds  is  plotted.  The  maximum  tension  is 
approximately  60,  000  lb  as  compared  with  38,  500  lb  for  the  case  of  6-ft 
waves  with  the  same  period. 

As  an  experiment  to  aid  in  understanding  the  effect  of  the  drag  caused 
by  the  presence  of  the  fluid  on  the  motion  of  the  cable,  one  case  was  carried 
out  with  zero  drag  (i.  e. ,  motion  in  vacuum).  A  very  interesting  motion 
pattern  was  obtained  which  appears  not  to  possess  a  periodic  character. 

This  solution  is  reproduced  in  Figure  8. 

The  successful  solution  of  this  problem,  as  well  as  a  number  of  others, 
involving  complex  nonlinear  systems  of  partial  differential  equations  by  the 
use  of  high-speed  calculators  and  finite  difference  methods  constitutes,  in 
the  opinion  of  the  authors,  a  major  advance  in  applied  mathematics. 
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Until  recently  it  has  been  considered  unt’easiole  to  obtain  numerical 
solutions  for  general  systems  of  partial  differential  equations  with  the 
exception  of  a  few  isolated  simple  types  of  equations  whose  solutions  are 
known  in  analytic  form.  However,  the  solution  of  engineering  problems, 
in  almost  every  major  field  of  science,  is  expressible  in  terms  of  systems 
of  partial  differential  equations.  Supersonic  and  subsonic  aerodynamics, 
nuclear  reactor  design  theory,  heat  flow,  propagation  of  electromagnetic 
and  acoustic  waves  are  but  a  few  areas  which  fall  in  this  category.  In  the 
past  engineers  have  largely  depended  on  experience  and  on  simplified 
linearized  models  of  the  phenomena  under  study.  In  the  future,  such 
simplified  theoretical  models  will  become  less  valid  -  as  speeds  under 
consideration  increase,  stresses  become  larger,  temperatures  higher. 

It  may  also  be  expected  that  experimentation  will  become  more  costly, 
more  time  consuming,  and,  at  times,  unfeasible.  It  is  fortuitous  that,  at 
the  same  time,  a  new  approach  appears  to  be  unfolding  for  the  solution  of 
many  difficult  engineering  problems  -  based  on  the  mathematical  represen¬ 
tation  of  the  phenomenon  and  the  numerical  solution  of  the  resulting 
unabridged  system  of  equations  by  the  use  of  high-speed  calculators  and 
finite  difference  methods. 

The  programming  of  the  various  phases  of  this  problem  was  carried 
out  by  Mr.  Thomas  Me  Fee,  of  the  Applied  Mathematics  Laboratory,  in  a 
most  effective  manner.  The  speed  and  accuracy  with  which  he  accomplished 
this  phase  of  the  solution  were  largely  responsible  for  the  success  in  meeting 
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the  required  time  schedules.  The  authors  would  also  like  to  express  their 
gratitude  to  Mr.  R.  T.  McGoldrick,  of  the  Structural  Mechanics  Laboratory, 
for  proposing  this  problem  and  for  a  number  of  helpful  discussions;  to 
Dr.  R.  Bart,  Structural  Mechanics  Laboratory,  for  a  numoer  of  ideas 
used  in  setting  up  the  numerical  procedure;  to  Dr.  E.  H.  Kennard,  David 
Taylor  Model  Basin,  and  Dr.  R.  M.  Langer,  Bureau  of  Ships,  for  helpful 
discussions  in  connection  with  the  definition  of  the  problem;  to  Dr.  Daniel  Shanks, 
Applied  Mathematics  Laooratory,  for  valuable  suggestions;  and  to 
Miss  Corinne  Lundgrep,  Applied  Mathematics  Laboratory,  for  assistance 
in  the  preparation  of  the  figures. 
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